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Abstract — Iterative image reconstruction (IIR) with sparsity- 
exploiting metliods, such as total variation (TV) minimization, 
investigated in compressive sensing (CS) claim potentially large 
reductions in sampling requirements. Quantifying this claim for 
computed tomography (CT) is non-trivial, because both full sam- 
pling in the discrete-to-discrete imaging model and the reduction 
in sampling admitted by sparsity-exploiting methods are ill- 
defined. The present article proposes definitions of full sampling 
by introducing four sufficient-sampling conditions (SSCs). The 
SSCs are based on the condition number of the system matrix 
of a linear imaging model and address invertibility and stability. 
In the example application of breast CT, the SSCs are used 
as reference points of full sampling for quantifying the under- 
sampling admitted by reconstruction through TV-minimization. 
In numerical simulations, factors affecting admissible undersam- 
pling are studied. Differences between few-view and few-detector 
bin reconstruction as well as a relation between object sparsity 
and admitted undersampling are quantified. 

Index Terms — Computed Tomography, Discrete-to-discrete 
Imaging Models, Sampling Conditions, Total Variation, Com- 
pressive Sensing, Breast CT 



I. Introduction 

RECENTLY, iterative image reconstruction (IIR) algo- 
rithms have been developed for X-ray tomography [I], 
Ja, |3i, lU, f5l, r6l, {7], f8l, f51, fTOl, fTTl based on the ideas 
discussed in the field of compressive sensing (CS) IIT2I . lfT3l . 
lfT4l . ifTSl . These algorithms promise accurate reconstruction 
from less data than is required by standard image reconstruc- 
tion methods. This is made possible by exploiting sparsity, 
i.e., few non-zeroes in the image or of some transform applied 
to the image. One can argue about whether these algorithms 
are truly novel or not: edge-preserving regularization and 
reconstruction based on the total variation (TV) semi-norm 
|[T6l, flTl, flS], I191 have a clear link to sparsity in the 
object gradient and have been considered before the advent 
of CS, and algorithms specifically for object sparsity have 
been developed for blood vessel imaging with contrast agents 
|^|. Nevertheless, the interest in CS has broadened the 
perspective on applying optimization-based methods for IIR 
algorithm development for computed tomography (CT), and 
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it has motivated development of efficient algorithms involving 
variants of the £i-norm fTTl, f2T1, flT^i, [23 1, 124}, ll26l . 

What is seldom discussed, however, is that the theoretical 
results from CS do not extend to the CT setting. CS only 
provides theoretical guarantees of accurate undersampled re- 
covery for certain classes of random measurement matrices 
lfT4l . not deterministic matrices such as CT system matrices. 
While the mentioned references demonstrate empirically that 
CS-inspired methods do indeed allow for undersampled CT 
reconstruction, there is a fundamental lack of understanding 
of why, and the conditions under which, this is the case. One 
problem in uncritically applying sparsity-exploiting methods 
to CT is that there is no quantitative notion of full sampling. 

Most IIR, including sparsity-exploiting, methods employ a 
discrete-to-discrete (DD) imaging modejj which requires that 
the object function be represented by a finite-sized expansion 
set and sampling specified over a finite set of transmission 
rays. This contrasts with most analysis of CT sampling, which 
is performed on a continuous-to-discrete (CD) model. For 
analyzing analytic algorithms such as filtered back-projection 
(FBP), a continuous-to-continuous (CC) model such as the X- 
ray or Radon transform is chosen, and discretization of the data 
space is considered, yielding results for the corresponding CD 
model. Analysis of the CD model is performed independent 
of object expansion. If the expansion set for the DD model is 
chosen to be point-like, e.g. pixels/voxels, there may be simi- 
larity between CD and DD models justifying some crossover 
of intuition on sampling, but in general sufficient-sampling 
conditions can be different for the two models. That a more 
fine-grained notion of sufficient sampling is needed for the 
DD model can be seen by considering the representation of 
the object function on a 128x128 versus a 1024x1024 pixel 
array. Clearly, the latter case requires more samples than the 
former, but sampling conditions derived from the CD model 
cannot make this distinction. Sufficient sampling for the DD 
model becomes even less intuitive for non-point-like expansion 
sets such as natural pixels, wavelets, or harmonic expansions. 
Yet, to quantify the the level of Mnt/ersampling admitted by 
a sparsity-exploiting IIR method, full sampling needs to be 
defined for the corresponding DD model, and to that end we 
introduce several sufficient-sampling conditions (SSCs). 

Specifically, in the present article, SSCs for the DD model 
are derived from the condition number of the corresponding 
system matrix. Multiple SSCs are defined to characterize both 

'See E2I Chapter 15] for an overview of different imaging models. 
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invertibility and stability of the system matrix. To perform 
the analysis, a class of system matrices is defined so that the 
system matrix depends on few parameters. The class is chosen 
so that it has wide enough applicability to cover thoroughly 
a configuration/expansion combination of interest, but not so 
wide as to make the analysis impractical. For the present study, 
we select a system matrix class for a 2D circular fan-beam 
geometry using a square-pixel array. The SSCs are chosen 
so that they provide a useful characterization of any system 
matrix class, but the particular values associated with the SSCs 
in this work apply only to the narrow system matrix class 
defined. While the article presents a strategy for defining full 
sampling, the analysis must be redone with any alteration to 
the system matrix class. 

After deriving the SSCs for the particular circular fan-beam 
CT system matrix class, we apply sparsity-exploiting IIR in the 
form of constrained TV-minimization. We consider the specific 
application of CT to breast imaging and use a realistic and 
challenging discrete phantom. We use the SSCs as reference 
points of full sampling for quantifying the undersampling 
admitted by each of the conducted reconstructions. Specifi- 
cally, we demonstrate significant differences in undersampling 
admitted for reconstruction from few views compared to few 
bins. We study how variations to the reconstruction optimiza- 
tion problem, to the image quality metric, to the discretization 
method for the system matrix, and to the sparsity of the 
phantom image affect the results. 

In Sec. [n] we describe the CT imaging model and present 
the particular system matrix class we employ for circular fan- 



beam CT reconstruction. In Sec. Ill we give a background 



on sparsity-exploiting methods. In Sec. IV the SSCs are 



presented and their application is illustrated for the 2D circular 
fan-beam case. Finally, Sec. |V] illustrates an example study 
on quantifying admissible undersampling by constrained TV- 
minimization employing the discussed SSCs. 

II. A CLASS OF SYSTEM MATRICES FOR THE 
DISCRETE-TO-DISCRETE IMAGING MODEL 

A. The X-ray transform 

Explicit image reconstruction algorithms such as FBP are 
based on inversion formulas for the CC cone-beam or X-ray 
transform model. 



dtf{s + t0), 



(1) 



where g, the line integral over the object function / from 
source location s in the direction 9, is considered data. Fan- 
beam FBP, for example, inverts this model for the case 
where the source location s varies continuously on a circular 
trajectory surrounding the subject, and at each s the ray- 
direction 9 is varied continuously through the object in the 
plane of the source trajectory. 

B. The discrete-to-discrete model 

For most IIR algorithms, the CC imaging model is dis- 
cretized by expanding the object function in a finite expansion 
set, for example, in pixels/voxels. Furthermore, the discrete 



digital sampling of the CT device is accounted for by directly 
using the sampled data without interpolation. The effect of 
both of these steps is to convert the imaging model to a 
discrete-to-discrete (DD) formulation. 



9 = Xf, 



(2) 



where g represents a finite set of ray-integration samples, / 
are coefficients of the object expansion, and X is the system 
matrix modeling ray integration. This DD imaging model is 
almost always solved implicitly, because the matrix X, even 
though sparse, is beyond large for CT applications: X is in 
the domain of a giga-matrix for 2D imaging and a tera-matrix 
for 3D imaging. 

A central point motivating the strategy of the present work 
is that the DD imaging model has narrower scope than the 
CD model, because it often derives from the CD model by 
expanding the continuous image domain with a finite set of 
functions. How the discretization of the CD model is done 
for CT to achieve the DD imaging model is not standardized. 
Many expansion elements have been used in CT studies; in 
addition to pixels, for example blobs |28|, wavelets ||29l , ll30l . 
and natural pixels lISTI . Il32l . Also, the matrix elements using 
only the pixel expansion set can be calculated in different ways 
that all tend toward the CC model in the limit of shrinking 
pixel size and detector-bin size. Different modeling choices 
will necessarily alter X. This tremendous variation in X 
means that it is important to fully specify X for each study, 
and it is important to re-characterize X for any change in the 
model. For example, changing pixel size can have large impact 
on the null space of the system matrix in the DD model. 

In order to describe precisely and provide a delimitation of 
the system matrices considered in the present work, we intro- 
duce the notion of a system matrix class. Any given system 
matrix depends on numerous model parameters determining 
the scanning geometry, sampling and discrete expansion set. 
A system matrix class consists of the system matrices arising 
from fixing a number of these parameters and leaving a subset 
of the parameters free. The system matrix class can then be 
studied by varying these free parameters. 

C. The system matrix class used in the present study 

In CT, projections are acquired from multiple source loca- 
tions which lie on a curve trajectory and the source location 
s(A) is specified by the scalar parameter A. The circular 
trajectory is the most common, and is what we use here, 

s( A) — Ro (cos A, sin A) , 

where Rq is the distance from the center-of-rotation to the X- 
ray source, and set to i?o = 40 cm in the present work. The 
detector bin locations are given by 

6(A, u) ~ (i?o — D){cos A, sin A) + u{— sin A, cos A), 

where D is the source-to-detector-center distance (D = 80 cm 
in the present work), and u specifies a position on the detector 
The ray direction for the detector-geometry independent data 
function is 

b(X, u) - s(A) 



9{X,u) 



|6(A,«)-,?(A)||2 
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The 27r arc is divided into A'views equally spaced angular 
intervals, so that the source parameters follow 



Xi = iAX, where 
AA = 27r/Arviews and i G [0, A^'views - !]• 



The detector is subdivided into M 



binsi 



(j + 0.5) At 



(3) 

(4) 



(5) 



where Dl is the detector length (Dl = 41.3 cm), Ujnin = 
~Dl/2, Au = Dl/Nuu,, and j G [0, A'bins - !]■ The detector 
length is determined by requiring it to detect all rays passing 
through the largest circle inscribed within the square N x N 
image array for which we use the side length 20 cm. We 
restrict the unknown pixel values to lie within this circular 
field-of-view (FOV), and the number of unknown pixel values 



pix 



(6) 



where the actual value, which has to be an integer, is given 
with each simulation below. 

Effectively, the dimensions of the projector X are M — 
-^views X -^bins rows (number of ray integrations) and iVpix 
columns (number of variable pixels). To obtain the individual 
matrix elements, the line-intersection method is employed, 
where Xm,n is the intersection length of the mth ray with 
the nth pixel. This description completely specifies the system 
matrix class for the present circular fan-beam CT study, and 
the free parameters of this class are N, A'views, and -/Vbins- 

III. CT IMAGE RECONSTRUCTION BY EXPLOITING 
GRADIENT-MAGNITUDE SPARSITY 

Reconstruction of objects from undersampled data within 
the DD imaging model corresponds to a measurement matrix 
X with fewer rows than columns. The infinitely many solu- 
tions are narrowed down by selecting the sparsest one, i.e., 
the one that has the fewest number of non-zeroes, either in 
the image itself or after some transform has been applied to 
it. Mathematically, the reconstruction can be written as the 
solution of the constrained optimization problem 



argmm 
/ 



*(/)||o such that Xf = g. 



(7) 



Here, ^I* is a sparsifying transform, for instance a discrete 
wavelet transform, and || • ||o is the ^o-"norm" (although it 
is in fact not a norm), which computes its argument vector's 
sparsity, that is, counts the number of non-zeroes. The equality 
constraint restricts image candidates to those agreeing exactly 
with the data. 

Central results in CS derive conditions on X drawn from 
certain random system matrix classes such that /* is exactly 
equal to the underlying unknown image that gave rise to the 
data g. Two key elements are sparsity of ^(/) and incoherence 
of X: exact recovery depends on the size of g being larger 
than some small factor of ||^'(/)||o IE], and the concept of 
incoherence is needed to ensure that the few measurements 
g available give meaningful information about the non-zero 



elements of ^'(/). Other important results in CS involve the 
relaxation of the non-convex £o-"iiorm" to the convex £i-norm. 



argmm 



|^'(/)||i such that Xf = g. 



(8) 



In contrast to Q, this convex problem is amenable to solu- 
tion by a variety of practical algorithms, although the large 
scale of CT matrices still presents a challenge for algorithm 
development. Another important contribution from CS is the 
derivation of conditions under which the solution to ([8]) is 
identical to the solution to (|7|, so that the sparsest solution 
can be found by solving (|8]l. 

For application to medical imaging, it was suggested in ifTSll 
that a potentially useful ^' would be to have 'i/ compute the 
discrete gradient magnitude of /, i.e., for the jth pixel 



*(/) . = ll^.-/l|2, 



(9) 



where Dj computes a finite-difference approximation of the 
gradient at each pixel j, and the 2-norm also acts pixel-wise 
on the differences. In CT, for example, the typical image 
consists of regions having an approximately constant gray- 
level value separated by sharp boundaries between various 
tissue types. The magnitude of the spatial gradient of such 
images is zero within constant regions and non-zero along 
edges, so the gradient magnitude image can be sparse. The 
£i-norm applied to the gradient magnitude image is known as 
the total variation (TV) semi-norm. 



|TV 



(10) 



(11) 



l*(/)l|l=Ell^^/1l2' 

i=i 

and the optimization problem of interest becomes 

/* — argmin ||/||Ty such that Xf — g. 
f 

However, the theoretical results from CS do not extend to 
the CT setting. Three properties that separate CT matrices 
from typical CS matrices are that CT matrices 

1) are structured and do not belong to random matrix 
classes for which CS results are proved |14|, 

2) can have rank smaller than the number of rows, which 
means that there exist vectors g in the data space that are 
inconsistent with X, and accordingly the linear imaging 
model (|2| has no solution, and 

3) may be numerically ill-conditioned in case of having 
more rows than columns (data set size is greater than 
the image representation). 

Nevertheless, it has been demonstrated empirically in exten- 
sive numerical studies with computer phantoms under ideal 
data conditions as well as with actual scanner data that highly 
accurate reconstructed images for "undersampled" projection 
data can be obtained from (|8]l and variants thereof. 

It is precisely this last phrase which is of interest in the 
present paper: what exactly does it mean to have "undersam- 
pled" data for CT? Under-sampled data implicitly relies on a 
certain level of sampling being sufficient — ^but no such precise 
concept exists for CT using the DD imaging model, to the 
best of our knowledge. Without a reference point for having 
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sufficient sampling it is difficult to quantify admissible levels 
of undersampling. In the present paper, we aim to provide 
this reference point. Specifically in the following section, we 
propose sufficient-sampling conditions (SSCs) to be computed 
for specific system matrix classes, and which serve as a 
reference for quantifying the admissible undersampling for re- 
construction with sparsity-exploiting methods. The application 
of the SSCs is demonstrated with numerical simulations of 
breast CT. 

IV. Sufficient-sampling conditions 

In considering sufficient sampling for circular fan-beam CT, 
the CC model is recast as a CD model by introducing a discrete 
sampling operator, usually taken to be evenly distributed delta 
functions, on the CT sinogram space. Making the assumption 
that the underlying sinogram function is band-limited, many 
useful and widely applicable results have been obtained, see 
for example Section 3.3 of ||33l and Refs. ||34|, ES, ll36l . 
Furthermore, for more advanced scanning geometries and 
sampling patterns there are available tools for analysis such 
as singular value decomposition (SVD), direct analysis of 
multi-dimensional aliasing, and the evaluation of the Fourier 
crosstalk matrix 1371 . Il38l . ||39l . These important results, 
however, do not apply directly to IIR, since for the DD model 
we need to take into account the finite image expansion set. 

We consider an empirical approach for characterizing suffi- 
cient sampling within a class of system matrices. The idea is 
to fix the image representation, which for the circular fan- 
beam system matrix class is the parameter TVpix, and then 
vary the sampling parameters A'views and A'bins to ensure 
accurate determination of the pixel values. This is done by 
establishing sufficient-sampling conditions (SSCs) based on 
matrix properties in the considered system matrix class. If 
the system matrix class is altered, the SSC-analysis must be 
redone. 

A. SSC definitions 

The SSCs, we propose, characterize invertibility and ill- 
conditioning of the system matrix class. In considering the DD 
imaging model (|2]), the data g are restricted to the range of 
X. This separates out the issue of model inconsistency which 
does not have direct bearing on sampling conditions. 

We define SSCl to be sampling such that X has at least 
as many rows as columns. If there are fewer rows than than 
columns, X necessarily has a non-trivial null space, and 
solutions of (|2]i will not be unique. Even if the number of 
rows is equal to or larger than the number of columns, there 
may still be a non-trivial null space, because the rows can 
be linearly dependent. In addition to SSCl, we define SSC2 
to mean that X has a null space consisting only of the zero 
image, or equivalently, that the smallest singular value a^in 
of X is non-zero. Existence of a unique solution to (|2| is 
ensured by SSC2. Both SSCl and SSC2 can be evaluated for 
any system matrix class. 

Neither of SSCl and SSC2 address numerical instability 
and to address that, we employ the condition number of X, 



the ratio of the largest and smallest singular values, 

nix) = O-max/o-min. (12) 

The condition number k can be as small as 1 and the larger k 
becomes, the more numerically unstable is solution of Xf = 
g. How to use k to define a SSC requires some discussion. 

Whereas sensing matrix classes studied in CS typically 
are well-conditioned — for instance the square discrete Fourier 
transform (DFT) matrix is orthogonal, thus having a condition 
number of 1 — the system matrices encountered in X-ray CT 
can have a relatively large condition number |33l, 1*401, which 
leads to numerical instability and thus large sensitivity to noisy 
measurements. Even if SSC2 holds, the condition number 
K is finite but may still be large, potentially allowing other 
images than the desired solution to be numerically close to 
satisfying Xf — g. If we fix the image representation, which 
for the present 2D circular fan-beam setup amounts to fixing 
iVpix, and increase the sampling, allowing iVyiews and A'bins to 
increase toward oo, the condition number will decrease toward 
a limiting condition number, 

kdc = lim k{X), (13) 

where the DC subscript refers to the fact that X is limiting 
to a discrete-to-continuous (DC) system matrix. The limiting 
condition number kjjc is the best-case k for a fixed image 
representation, but kdc may still be larger than 1. 

For actual CT scanners, it is not practical to allow A'views 
and A'bins to increase without bound, and empirical experience 
shows diminishing improvements in doing this. To balance 
the impracticality of going to continuous sampling on the one 
hand against the need to optimize numerical stability on the 
other, we introduce SSC3 to mean that the condition number 
of X satisfies 

k{X)/kdC < ?'samp, (14) 

where rsamp is a finite ratio parameter greater than 1. The 
smaller the choice of rsamp, the closer X is to the DC limit. 
This SSC can also be generally applied to other system matrix 
classes, but the appropriate parameter setting of rsamp will be 
specific to a particular class. 

Finally, we introduce SSC4 specifically for the present 2D 
circular fan-beam system matrix class. This SSC is taken to 
mean 2N samples in both the view and bin directions, i.e., 
-Aviews = Abins = 2A'. This SSC is simple to evaluate, and 
we will demonstrate empirically that it is a useful condition, 
which acts as a good approximation for attaining SSC3 with 
''samp = 1-5. This SSC is specific to the system matrix class 
investigated here. Even slightly different system matrix classes 
might not allow for the same SSC4 definition. 

Our strategy is similar to analysis presented in early works 
on CT, such as in |40|, but the point here is not novelty of 
the analysis; rather we need to establish a reference point by 
which to evaluate the sampling reduction admitted by sparsity- 
exploiting methods. 

In what follows, the proposed SSCs are examined for the 2D 
circular fan-beam system matrix class. First, small systems are 
considered, where X can be explicitly computed and analyzed 
so that the full set of singular values of X is attainable. 
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Second, we argue that our conclusions generalize to larger, 
more realistic systems, where X is impractical to store in 
computer memory and it is only feasible to compute the 
smallest and largest singular values. 

B. SSCs for small systems 

We consider a small N — 32 image array with iVpix = 
812, and generate system matrices X for different numbers of 
views, iVviews € [8,128], and detector bins, iVbins G [8,128]. 
The condition number k{X) is computed through direct SVD 
of X for all values of A'views and A^bins within the specified 
parameter ranges, and k{X) as a function of A'views (for fixed 
-^bins = 64) and A'bins (for fixed A'views — 64) is shown in Fig. 
[T] The plots show three phases: the left-most part, where the 
condition number is infinite; the middle, where the condition 
number becomes finite and decays slowly; and the right-most 
part, where it remains relatively stable. The positions of the 
different SSCs are shown at the top and bottom, and they serve 
as transition points between the three phases. 

For varying the number of views, SSCl occurs at A'views — 
13, where X is of size 832 by 812. In fact, also SSC2 occurs 
here, since k has become finite. For varying the number of 
bins, SSCl occurs at the same place, but SSC2 needs 14 bins. 
In general, we have no way to determine whether SSCl and 
SSC2 occur in the same position for the whole system matrix 
class, which makes SSCl less reliable as a general reference of 
full sampling. On the other hand, SSC2 is a reliable reference 
point for full sampling, however, SSC2 requires more work to 
determine, because its location can change with a change of 
system matrix class. 

After passing SSC2, the condition number k decreases. For 
larger A'views and Abins, the decay becomes slower, and we pick 
''samp = 1-5 as a trade-off between a sufficiently small condi- 
tion number and a finite number of views. As an approximation 
of kdc we take the value of k at A'^^iews = A'bins = 4A' = 128, 
yielding kdc = 9.17. Then SSC3 occurs at A'views — 68 and 
at Abins = 68, which suggests a symmetry in A'views and Abins- 
On the other hand, the decrease in k during the middle part is 
not symmetric in the parameters A'views and Abins; the decrease 
in K with A'views is gradual while that of Abins is step-like at 
Abins = 48. Nevertheless, at the position of SSC3, there is 
only small further reductions in k to be gained by going to 
larger A'views and Abins- The simpler condition SSC4 occurs at 
^views — 64 and A'bins = 64, and it approximates SSC3 with 
J-samp = 1-5 closely. 

C. Altering the system matrix class 

Altering the system matrix class will in general alter 
the SSCs. To demonstrate this effect, we replace the line- 
intersection based system matrix class by ray-tracing, using 
nearest-neighbor interpolation at the mid-line of each pixel 
row. The experiment is repeated and the obtained condition 
numbers are shown in Fig. [2] along with the ones based on 
line-intersection, for comparison. The shown SSCs are for 
ray-tracing. While the same overall trends are seen, there 
are some significant differences. First, for fixed A'views = 64 
we need Abins — 16 to obtain SSC2, compared to 14 for 



line-intersection. This firmly establishes that SSCl does not 
imply SSC2, and that the precise position of SSC2 cannot 
be inferred from knowing SSC2 of a similar system matrix. 
Second, the ray-tracing condition numbers are smaller than the 
same for line-intersection, for instance, koc — 7.23 is 20% 
lower relative to the line-intersection version of X. That the 
ray-tracing condition numbers are lower does not necessarily 
mean that this method is "better" than the line-intersection 
method for real-world applications, because the other side of 
the story is model error, which is not considered here. Finally, 
the positions of SSC3 are different, for fixed Abins ~ 64 
coinciding with SSC4, while for fixed A'views — 64 occurring 
at Abins = 56. Still, for larger A'views and Abins there are 
only small further reductions in k, to be gained, and SSC4, 
at Tsanip = 1.45, approximates SSC3 closely. 

One could imagine that other system matrix classes such 
as employing area-weighted integration instead of the linear 
integration or different basis functions could alter the condition 
numbers of X even more substantially. We do not include 
results for more system matrix classes, as our goal here 
is not to provide a comprehensive comparison between all 
conceivable classes, but merely to stress that different classes 
can have different SSCs, and to propose carrying out the same 
study for gaining insight in the particular system matrix class 
at hand. 

D. SSCs for larger systems 

The results shown in Fig. [T] and Fig. |2] give a sense 
about various sampling combinations, but the system size is 
unrealistically small. In this section we aim to extend the 
results to larger systems. For large X, it is not practical to 
compute the direct SVD for evaluating k. Instead, we seek 
only to obtain (Tniin and dniax, which can be accomplished 
through the power and inverse power methods |41J. 

For characterizing the present circular, fan-beam system 
matrix class, k,{X) is computed for larger image arrays with 
sizes A'' = 32, 64, 128, 256. We focus on sampling conditions 
where A'bins = 2A'" and report the data sampling in views, 
-^viewsi as multiples of A', ranging from 1.0 to 4.0. The left 
plot in Fig. |3] shows the condition number as a function of 
A' for each sampling size on a double logarithmic scale. A 
clear linear trend is seen in all cases and the best Unear fits 
and their slopes are also shown, in all cases very close to 
0.50, and we conclude that k scales with \/N . For increasing 
A'views^ the condition numbers tend towards the bottom line, 
note in particular that not much difference is seen between 
-^views — 3N and A'views = 4Af indicating that the limiting k is 
approached. We conclude that kdc also scales with ^/N. As a 
result, it can expected that at SSC4, i.e., A^views = Abins = 2A^, 
k{X)/kjjc ~ 1-5, which was the case at A' = 32. Hence, 
SSC4 will continue to approximate SSC3 closely, when the 
image size is increased. To further support this conclusion we 
show in the right plot of Fig. |3]the ratio rsamp = kdc 
as function of the number of views (normalized by A') for 
each N . The r^^^^ values are almost identical for all A^ and 
intersect the line rsamp = 1-5 very close to A'views — 2A', which 
is precisely SSC4. 
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Fig. 3. Left: condition numbers as function of image size. Each symbol represents different number of views ranging from to 4A'^. Circular image arrays 
are used with sizes given by N = 32, 64, 128, 256. The number of bins is fixed at 2N. With each number of views is also shown the best linear fit and its 
slope is given. In all cases the condition number scales with \/7V. Right: same condition numbers normalized by the respective kdc at each N and plotted 
as function of view number normalized by image size A'^. The full line is the point-wise mean and the dashed line is the position of SSC3 with rsamp = 1-5. 
Independently of N, SSC3 with rsan,p = 1.5 occurs very close to A'^vicws = 2Af. 



E. Summary of SSCs 

The conditions SSCl and SSC2 are useful reference points 
for invertibility of X and can be computed for any system 
matrix class. The size of the gap between SSCl, X being 
square, and SSC2, X having an empty null space, is governed 
by inherent linear dependence of the rows of the system 
matrix. Because the results show little difference between 
SSCl and SSC2 for the present system matrix class and SSCl 
is easier to compute, we use only SSCl in the simulation 
studies in the following section. 

For system matrix classes representing CT imaging, stability 
of the system matrix plays an important role, and accordingly 
we have introduced SSC3 which also can be computed for any 
system matrix class. For the present, circular fan-beam system 
matrix class, SSCS at rsamp = 1.5 is a useful operating point, 
and this level of sampling is well approximated by the simple 
rule, SSC4, where iVyiews = ^'^bins = 2iV. We point out that for 
other system matrix classes, even those representing circular 
fan-beam CT, other operating points for SSC3 may be more 
appropriate and empirical studies must be performed to see if 
a simple condition, such as SSC4, can approximate accurately 
SSC3. 

In the remaining part of the paper we demonstrate how 
we can use the SSCs as a reference for stating admissible 
undersampUng factors in sparsity-exploiting reconstruction. 

V. Numerical experiments with SSCs and 

SPARSITY-EXPLOITING UNDERSAMPLED RECONSTRUCTION 

In this section we investigate sparsity-exploiting IIR in nu- 
merical simulation studies. Our goal is to numerically demon- 
strate and quantify the undersampling admitted by sparsity- 
exploiting IIR, i.e., at which an accurate reconstruction is 
obtained. We use numerical simulation, because we are un- 
aware of any theoretical results establishing undersampling 
guarantees for the present system matrix class. We focus 



here on exploiting gradient-magnitude sparsity by use of 
constrained TV-minimization. 

Three important factors differentiate the present stud- 
ies from previous simulation work with constrained TV- 
minimization: 

1) use of phantoms with realistic complexity, 

2) numerically accurate solution to the constrained TV- 
minimization problem, and 

3) quantitative references for full sampling — the central 
topic of the paper 

For each factor, we we briefly discuss the significance. 

Much simulation work on constrained TV-minimization has 
used regular, piece-wise constant phantoms, such as the Shepp- 
Logan phantom, to demonstrate the promise of the technique. 
For that purpose, such unrealistically simple phantoms were 
fine, and simulations were generally followed up by demon- 
stration with actual CT projection data. For the present purpose 
of quantifying admissible undersampling, we need phantoms 
with similar complexity as would be encountered in CT 
applications, and as an example we focus on breast CT. The 
standard measure of complexity employed in CS is the image 
sparsity, i.e., the number of non-zeroes in the image, or in 
the case of TV-minimization the number of non-zeroes in the 
gradient-magnitude image. Accordingly, we choose a digital 
phantom with realistic gradient-magnitude sparsity modeling 
breast anatomic structure ll42l . 

The accuracy requirement on the solver of constrained TV- 
minimization for the present study is extremely high. The 



optimization problems in Sec. V-B below, are solved to high 



accuracy, which has been made possible only recently for 
large-scale CT problems involving the TV-semi-norm through 
development of advanced first-order methods 1221 . Il23l . Il26l . 
This level of accuracy is necessary, because empirical image 
error results obtained by sweeping parameters of the system 
matrix class will be used for determining whether a numeri- 
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cally computed solution is close to the original. High-accuracy 
solutions remove any doubt about whether the resulting im- 
ages, and the corresponding quantitative measures, depend on 
the algorithm used to solve constrained TV-minimization. 

The SSCs defined above provide reference points useful 
for interpreting the empirical results of this section and 
help to quantify undersampling admitted by constrained TV- 
minimization. 

A. Breast CT background 

Breast CT B3| . Il44l is being considered as a possible 
screening or diagnostic tool for breast cancer. The system 
requirements are challenging from an engineering standpoint, 
because this type of CT must operate with a total exposure 
similar to two full-field digital mammograms (FFDM). FFDMs 
for a screening exam entail two X-ray projections, while 
breast CT acquires on the order of 500 X-ray projections. The 
exposure previously used for only two views is now divided 
up among 250 times more projections. Accordingly, sparsity- 
exploiting IIR algorithms for CT may have an impact on the 
breast CT application. The potential to reconstruct volumes 
from fewer views than a typical CT scan might allow an 
increased exposure per view. 

For the present study, we employ the breast phantom 
originally described in ll42l and displayed in Fig.|4] It consists 
of A^pix = 51,468 pixels within the circular image region, 
contained in a 256 x 256 array. The breast phantom has a 
small region of interest (ROI) containing 5 tiny ellipses which 
model microcalcifications. The gray values range from 1.0 
to 2.3. The modeled tissues and corresponding gray values 
are fat at 1.0, fibroglandular tissue at 1.10, skin at 1.15, and 
microcalcifications ranging from 1.9 to 2.3. The sparsity in 
the gradient magnitude image is 10, 019, or roughly one fifth 
of iVpix. Because we are investigating the utility of gradient- 
magnitude sparsity-exploiting algorithms, it is important that 
the test phantom have a realistic sparsity level relative to the 
actual application. 

B. Simulation optimization problems and algorithms 

Our goal is to evaluate quantitatively what level of un- 
dersampling reconstruction through ([TTJ allows. Similar to 
the analysis of the linear imaging model (|2|, only data g 
in the range of X is considered. Although not a realistic 
assumption for actual CT data, this "inverse crime" scenario 
P31 is appropriate for obtaining a reference of the underlying 
admissible undersampling. For the numerical studies, we solve 
a relaxed form of ( [TT| , where the data equality constraint is 
replaced by an inequality allowing for a small deviation e from 
data as measured by the distance D between the data g and 
the projection Xf of some image /, 

DiXf,g)<e, (15) 

where 

D{Xf,gf = \\Xf~g\\l 

* views bins 

Scaling the data error D with iVviews and A'bins is done to enable 
comparison across images reconstructed from different view 



and detector-bin numbers. The constrained TV-minimization 
problem is 

VTV: f =argmin||/||Tv (16) 
/ 

subject to D{Xf,g)<e. 

Accurate solution of ([T6| is non-trivial; although the objec- 
tive is convex, it is not quadratic. The algorithm employed here 
solves its Lagrangian using an accelerated first-order method, 
using only the objective and its gradient, and is explained in 
detail in [22 1. An important technical detail for this algorithm 
is that it requires that the image TV-term be differentiable. 
For the alg orithm imple mentation we use a smoothed TV- 
term, 'J2j y \\Djf\\2 + V' with a small smoothing parameter, 
T] = 10^ One convergence check on the algorithm is 
performed by evaluating 

(v/i^(/))-(v/W,g)) 

cosa=^ ^ — ^ -, (17) 

|V^~i?(/)||V^-D(X/,5)| 

where R{f) denotes a generic regularization term, and for 
constrained TV-minimization R{f) — W/Wtv- The conditions 
for convergence, derived in Ref. |4|, are that the gradients 
of the data-error and regularization terms are back-to-back, 
cos a = —1, and D{Xf,g) — e. The latter condition assumes 
that the data-error constraint is active, which is the case for 
all the simulations performed here. For the present results, 
iteration is terminated when both 

cos a < -0.9999 (18) 
|i^(X/,g)-e|/e< 0.001 

are satisfied. 

C. Admitted undersampling by i2-TV 

We are interested in two separate notions of accurate re- 
construction: exact reconstruction and stable reconstruction. 
By the former, we mean that the reconstructed image is 
identical to the original. Exact reconstruction is only possible 
when e = 0, because the regularizing effect of a non-zero e 
introduces a bias relative to the original image. Having e = 
is only relevant for noise-free data, which means that stability 
is not an issue. Since SSC2 ensures a unique solution to 
it can be used as a full sampling reference point for exact 
reconstruction. In practice, for the considered system matrix 
class, we found little difference in locations of SSCl and 
SSC2, and we will therefore use SSCl as a surrogate for SSC2. 

Stable reconstruction is the corresponding concept for fixed, 
non-zero e, where we cannot hope for exact reconstruction. 
Instead, we are interested in the degree of sampling at which 
further increase in sampling leads to no further improvement 
in the reconstruction. This point of stable reconstruction can 
be compared to SSC3, since that is the point where no further 
improvement of the system matrix condition number occurs. 
Since SSC4 was seen to approximate SSC3 for the present 
system matrix class and is simple to determine, it could also 
be considered to use SSC4 instead for the reference point. 
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Fig. 4. Left: 256 X 256 pixelized breast CT phantom used in the present study. Middle: same with ROI around microcalcifications shown magnified as inset. 
Right: the gradient magnitude image, which has a sparsity of 10, 019 non-zero pixel values. 
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In the simulations, we take the image array to be the same 
as that of the breast phantom TV = 256. The parameters of the 
circular, fan-beam system matrix class are varied in a fashion 
parallel to the condition number plots of Fig. [T] first, A'bins is 
fixed at 2N and A'view.s is varied in the range [32,512]; and 
second, A'^^iews is fixed at 27V and iVbins is varied in the range 
[32,512]. 

Computing SSCl and SSC4 is straightforward and they 
occur at A'views — 101 and iVviews — 512, respectively. Com- 
putation of SSC3 was performed using the procedure out- 
lined in Sec. |IV-D| For the fixed-bin case, SSC3 occurs at 
-^views = 492, with Tsamp = 1.51 and for the fixed-view case 
at iVbins = 456. 

In order to assess the undersampling with respect to ex- 
act reconstruction admitted by exploiting image gradient- 
magnitude sparsity, we need access to the solution of £2- 
TV for a data-equality constraint, e = 0. To our knowledge, 
solving i!2-TV for e = accurately with the present system 
size is currently impractical. Instead we solve ^2-TV for 
e = 10^"^, 10^^, and 10^^ to study the reconstruction error 



as e approaches zero. As an error measure, we use the root- 
mean-square-error (RMSE), 



'^ = iir-/oii2/. 



(19) 



where /* is the solution to ^2-TV and /□ is the original 
phantom. 

The computed RMSEs for the results from ^2-TV are 
displayed in Fig. |5] As in Sec. |IV] we show SSC-locations at 
the top and bottom. The horizontal line shows the minimum 
gray level contrast, 0.05, in the test phantom and provides 
a reference for the RMSE. An image RMSE much less 
than the minimum gray-level contrast is an indicator that the 
reconstructed image is visually close to the original phantom 
(bailing pathological distributions of the image error). 

For the plots of S versus A^views, we note a steep drop in 6 
as A'view.s increases past 40 views and the drop is increasingly 
rapid as e decreases. Based on these curves we extrapolate 



that exact reconstruction would be attained for N^i 



50 



at e 0. Because SSCl occurs at A'views = 101, we note an 
admitted undersampling with respect to exact reconstruction 
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of a factor of 2 for the present simulation. Note that use of 
SSCl leads to a conservative estimate, because SSC2 can only 
be larger than SSCl. 

For the plots of S versus iVbins, the image RMSE curves drop 
much more gradually at each of the e's investigated. Based on 
these curves it is only clear that S is tending to zero at A^bins = 
190 as function of e. Comparing to SSCl at A'bins = 101, we 
do not observe any level of admitted undersampling in the bin- 
direction with respect to exact reconstruction. Extending the 
range of e to smaller values may yield a different conclusion. 

This difference reflects an asymmetry in sampling of the 
two parameters of X. We note that the asymmetry was also 
observed in the condition number dependence on A'views and 
^bins for the iV = 32 simulations in Sec. IV-B| For the present 
N — 256 simulations, a relatively large k{X) = 3.2 • 10^ 
is seen for iVbins — 128 and A^views = 512, compared to 
k{X) = 1.5 • 10^ for iVbins = 512 and Ny,^^, = 128. The 
results demonstrate a larger potential for successful TV-based 
reconstruction from few views compared to using few bins. 

Regarding stable reconstruction, we note that the curves in 
Fig. |5] all exhibit a plateau, where 6 levels off with increasing 
A'view.s or iVbins, meaning that no gain in image RMSE is 
achieved by increasing sampling. Thus, the left-most point of 
these plateaus is the point of stable reconstruction. For the 
plot varying iVviews, we see that stable reconstruction begins 
at A^views ~ 80, which is a factor of 6 fewer than SSC3. 
For the plot varying iVbins, the stable reconstruction begins 
at A'bins ~ 200, a factor of approximately 2 fewer than SSC3. 

These results show quantitatively that significant undersam- 
pling in iVviews, particularly with respect to stable reconstruc- 
tion, is admitted for ^2-TV. This conclusion is achieved with 
a phantom modeling a realistic level of gradient-magnitude 
sparsity. We do point out, however, that these empirical 
results only apply to the presented simulation. To support the 
present conclusion for admitted undersampling, we vary in the 
following sections different aspects of the £2 -TV study. 



D. Altering the optimization problem 

To support the use of the gradient-magnitude sparsity 
exploiting i'2-TV for admitting undersampling, we compare 
results with two other optimization problems. 



and 



£2-magnitude: 
subject to 

^2-roughness: 
subject to 



/* = argmin ||/|| 
/ 

D{Xf,g)<e, 



f* = argmin II V/l 
/ 

D{X,fJ)<e, 



(20) 



(21) 



where V represents a numerical gradient operation and is com- 
puted by forward finite-differencing. The Lagrangian form of 
these optimizations are two forms of Tikhonov regularization 
commonly used for IIR. 

The solutions to ^2-magnitude and ^2-roughness are ob- 
tained with linear conjugate gradients (CG) applied to the 
Lagrangians of these problems with the multiplier being 
adjusted until the data-error constraint holds with equality. 
The convergence criteria are the same as what is specified 
in Eq. (I81 except that R{f) = II/II2 for £2-magnitude and 



R{f) = l|V/||i for ^2-roughness. 

We focus on the e = 10^^ case and plot image RMSEs for 
£2 -magnitude, £2-roughness and ^2-TV as function of A'views 
for fixed A'bins = 512 in Fig. [6] We note little difference 
between results from ^2-inagnitude and £2-roughness but a 
large gap between these results and those of ^2-TV. The 
optimization problems ^2-TV and €2-roughness differ only 
on the norm of the image-gradient in the regularization term, 
while ^2 -roughness and ^2 -magnitude differ by the presence 
of the gradient. It is clear from Fig. [6] that for this simulation, 
the £i-norm has the greater impact. 

For large A^^^iews, the £2 -TV RMSE is actually slightly 
inferior to that of ^2-magnitude. The reason is the regular- 
izing effect of having a non-zero e, which causes a small 
bias of the solutions compared to the original image. The 
relative size of the biases are not known in advance. We 
conclude that the £2-TV solution is not to prefer over the £2- 
magnitude and ^2-roughness solutions when A^viewii approaches 
the SSC3 (Tsamp — 1-5). Nevertheless, there is a certain 
"sampling window", for the present phantom, approximately 
for A^views G [50,256], where the constrained TV- minimum 
solution is superior to Tikhonov regularization in terms of 
RMSE. 

In Fig. [tJ we overlay the results of ^2-magnitude onto the 
results of ^2-TV from Fig. |5] in order to investigate possible 
undersampling admitted by ^2-magnitude. The results of £2- 
roughness are not shown because they are similar to those 
of ^2 -magnitude and to prevent clutter in the figure. Going 
from left to right, both plots show a gradual decrease of 6 
for ^2 -magnitude as A'views and Abins increase with 6 leveling 
off at A'view.s ~ 300 and Abins ~ 400. For the investigated 
range of e, £2 -magnitude does not admit any undersam- 
pling with respect to exact reconstruction, but does show a 
marginal undersampling with respect to stable reconstruction 
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as the corresponding (5-curves reach the plateau before SSC3 
and SSC4. In summary, the undersampHng admitted by £2- 
magnitude is substantially less that that admitted by 
for this simulation; particularly in considering view number 
undersampHng with respect to stability. 



E. Altering the image evaluation metric 

Conclusions based on evaluating reconstructed images with 
a single summarizing metric, such as the RMSE, can be 
misleading. While our aim is not a fully realistic image 
evaluation, we want to show how the results can potentially 
change with a change of metric. For example, with the task 
of microcalcification detection in mind, one might consider 
the RMSE of only the ROI of the microcalcifications dis- 
played in Fig. |4] This RMSE is denoted 5-g.oi- We show 
the corresponding plot for (5roi in Fig- [8] While there are 
numerical differences between the 5 and (5roi plots, the trends 
are similar giving us further confidence that the RMSE of 



the entire image, 5, can be used for investigating admitted 
undersampHng. 

Another way to evaluate images is by visual comparison. 
The reconstructed images in Fig. [9] are shown for a range in 
A'views showing the transition to accurate image reconstruction 
by £2-TV. That the €2-magnitude and ^2-roughness show 
strong artifacts for this range is expected as their coiTe- 
sponding image RMSEs are at the level of the minimum 
phantom contrast level. Interestingly, the microcalcifications 
can be identified and well-characterized in all reconstructions, 
although more clearly with more views. It may be argued that, 
from a utility point of view, that 32 views would suffice if we 
are solely interested in the microcalcifications and disregard 
the prominent artifacts of the background image. The ROI of 
the ^2 -TV reconstructed images are shown with a narrow gray 
scale window in the bottom row to reveal the high level of 
accuracy at A'views = 50. We emphasize here that our goal 
is not to go into a discussion about different artifacts but 
simply support our conclusions on undersampHng from Sec. 



V-C by illustrating the behavior in the transition region around 

A'view.s ~ 50. 

F. Altering the system matrix class 

To illustrate the change in results due to change in the 
system matrix class, the RMSE, 5, is again computed as 
function of iVvjews for Nh\n^ — 2N 
10-'^ and 10 

ray-tracing with nearest-neighbor interpolation at the mid 



views for iVbins ^ 2N = 512 and for e = 10"^, 
but using a system matrix set up through 



line of each pixel row, as in Sec. I V-C Results are plotted 
in Fig. 10 The overall trends are similar to those for line- 
intersection, however SSC3 occurs already at A'views — 408, 
with Tsamp = 1.49. By closer comparison with the line- 
intersection results, it is seen that the nearest-neighbor RMSEs 
are smaller than the line-intersection RMSE at the same A'views- 
This example serves to illustrate that the SSCs will change 
when the system matrix class is altered. 

While the present alteration is relatively minor, it is enough 
that the approximation SSC4 of to SSC3 is worse, and larger 
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Fig. 9. 1st, 2nd, 3rd and 4th columns show reconstructions from 32, 40, 48, 64 view data. The 1st, 2nd and 3rd rows show ^2-magnitude, £2-roughness and 
£2-TV images with e = 10""". The gray scale window is [0.95, 1.20] for the complete image, and [0.9, 1.8] for the ROI blow-ups. The 4th row shows the 
^2-TV ROIs enlarged in an extremely narrow gray scale, [1.09, 1.11], in order to scrutinize the transition to sufficient sampling based on the object sparsity. 
These images show that 50 views is sufficient for this system and object. 



differences can be expected with more radical changes such 
as the use of non-point-like image expansion functions. 



In terms of admissible undersampling for ^2 -TV with the 
altered system matrix class, we see very similar undersampling 
factors for both exact and stable reconstruction as for the line- 
intersection class. 



G. Altering the phantom sparsity 

The breast phantom study is repeated employing a variation 
of the FORBILD head phantom [46 1 which is highly sparse 
in the gradient-magnitude image. The present version of the 
phantom, which is seen in Fig. 



11 does not have the ear 



objects of the original phantom, and the contrast levels have 
been increased so that the minimum gray-level contrast is 
the same as for the breast phantom. The gradient magnitude 
sparsity is 2.492, or approximately a quarter of the breast 
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SSC3 SSC4 ensembles of phantoms in order to verify that admitted under- 

sampHng for constrained TV- minimization depends primarily 
on the gradient-magnitude image sparsity. 

We also note that while we may have exact reconstruction 
of the head phantom at A'yiews = 12 and the reconstructed 
image at e = 10^^ appears very accurate in the [0.9, 1.1] gray 
scale window, it is in fact not an exact reconstruction. By 
nan^owing the gray scale to [0.99, 1.01], also shown in Fig. 
[TT] prominent artifacts become visible. This underlines that 
exact reconstruction is not the relevant notion for a fixed, non- 
zero e. Instead, stable reconstruction, at iVviews = 64, yields 
an accurate reconstruction, and for the present case already at 
512 ^views = 32 (also shown in Fig. [TT| i the artifacts are reduced 
to a negligible level in the [0.99, 1.01] gray scale window. 



Fig. 10. Image RMSE, S, for £2 -magnitude and £2-TV reconstructions, as 
in Fig. |6] for varying e, computing X witli ray-tracing and nearest-neighbor 
interpolation. A^bins is fixed at 512. The labels SSCl, SSC3 and SSC4 are the 
sufficient-sampling conditions discussed in Sec. |IV-A| 



phantom. In Fig. [TT] the obtained 6 for ^2 -magnitude, £2- 
roughness and £2 -TV are shown as function of iVyiews with 
iVbins = 2N = 512 and e = 10"^ 

The ^2-magnitude and ^2-roughness curves are almost iden- 
tical to those of the breast phantom. The SSCs are unchanged, 
because the same system matrix class is used. That two 
so different looking phantoms show such similar (5-behavior 
suggests that the reconstruction quality of ^2-magnitude and 
^2-roughness depend only weakly on the particular phantom. 

For £2 -TV, on the other hand, the step-like transition occurs 
already at iVyiew., = 12, for which the reconstruction is shown 
We expect that this phantom would be recovered 
M/iews = 12 in the limit e 0, leading to an 



in Fig 11 



exactly at 

admitted undersampling with respect to exact reconstruction 
of a factor of approximately 8. Stable reconstruction occurs 
at iVviews ~ 64, i.e., an undersampling also of a factor 8 with 
respect to stability. 

Interestingly, the exact reconstruction result hints at the 
existence of a simple relation between sparsity and admitted 
undersampling. Compared to the breast phantom, there is a 
gain in undersampling by 8/2 4. In comparison, we note 
that the change in gradient magnitude sparsity relative to the 
breast phantom is 10,019/2,492 « 4. That is, reducing the 
sparsity by a certain factor leads to an improvement in the 
admitted undersampling by the same factor. This result, if 
shown to hold, could be important for practical application 
of CS-inspired sparsity-exploiting methods, since it provides 
quantitative insight into how many views would suffice for 
reconstructing images of given sparsity. Another conclusion 
that can be drawn is that simulations with images of too low 
sparsity compared to a realistic level in the imaging scenario 
of interest are bound to yield over-optimistic promises of 
undersampling potential. This could have been anticipated but 
the result establishes this expectation quantitatively. 

We caution, however, that the result is based on only two 
phantoms and further study is required. For more robust 
conclusions, the present studies need to be performed on 



VI. Summary 

We argue that a quantitative notion of a sufficient-sampling 
condition (SSC) for X-ray CT using the DD imaging model 
is necessary in order to provide a reference for evaluating 
the undersampling potential of sparsity-exploiting methods. 
We propose and apply four different SSCs to a class of 
system matrices describing circular, fan-beam CT with a pixel 
expansion. While SSCl and SSC2 characterize invertibility of 
the system matrix, SSC3 characterizes numerical stability for 
inversion of the system matrix. A simple-to-compute SSC4 is 
seen to approximate SSC3 closely for the circular, fan-beam 
full angular range CT geometry. 

We employ the SSCs as reference points of full sampling 
to quantify undersampling admitted by reconstruction through 
TV-minimization on a breast CT simulation. Relative to SSCl, 
we observe some undersampling potential of TV-minimization 
for exact reconstruction and large undersampling relative to 
SSC3 and SSC4 for stable reconstruction. We find few- 
view reconstruction to admit larger undersampling than few- 
detector bin reconstruction, and we show evidence of a simple 
quantitative relation between image sparsity and the admitted 
undersampling. 

More generally, the proposed SSCs can help to engineer 
and understand other sparsity-exploiting IIR algorithms by 
providing full sampling reference points for the system matrix 
class associated with the imaging application of interest in 
order to quantify the admissible undersampling. This analysis 
can guide decisions on alternative optimization problems, ob- 
ject representations, sampling configurations, and integration 
models. 
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Fig. 11. Top: Image RMSE, S, for ^2-magmtude, ^2-roughness and ^2-TV reconstructions, as in Fig.|6] except that the data are generated from the head 
phantom shown on the right. The labels SSCl, SSC3 and SSC4 are the sufficient-sampling conditions discussed in Sec. |IV-A| Bottom: images reconstructed 
for A'^vicws = 12 shown in a gray scale window of [0.9, 1.1] on the left and a narrow gray scale window of [0.99, 1.01] in the middle. On the right is the 
reconstructed image for A^vicws = 32 in the naiTow gray scale window [0.99, 1.01]. At Afyicws = 12 the RMSE is 0.005 resulting in visible artifacts for the 
image shown in the 1% gray scale window, while the RMSE is a factor of 10 less for A^vicws = 32. 
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